home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / T-test_IndSamples.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1998-08-15  |  10.4 KB  |  547 lines

  1. /* T-Test 2 Means for Independent Samples */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Math Library needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30) 
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19.  /* add to library list */
  20. signal off syntax
  21.  
  22. /* Start Main Routine */
  23. if loadflag=1 then 'Load()'
  24. 'ActivateWindow()'
  25. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  26. colon=pos(":",range)
  27. if colon=0 then do
  28.     'Message "Please select a range before executing this script"'
  29.     'DEFPUBSCREEN("Workbench")'
  30.     exit
  31. end
  32.  
  33. /* Find cell references and cell, column numbers */
  34. start_cell=substr(range,1,colon-1)
  35. end_cell=substr(range,colon+1)
  36. start_row=cellrow(start_cell)
  37. end_row=cellrow(end_cell)
  38. start_col=cellcol(start_cell)
  39. end_col=cellcol(end_cell)
  40. NRows=end_row-start_row+1
  41. NCols=end_col-start_col+1
  42. if NCols >2 then do
  43.     'Message "Only Two Samples Allowed"'
  44.     'DEFPUBSCREEN("Workbench")'
  45.     exit
  46. end
  47. /* Get cell reference for output range */
  48. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  49. if out_cell="" then do
  50.     'DEFPUBSCREEN("Workbench")'
  51.     exit
  52. end
  53. if length(out_cell)<2 | datatype(left(out_cell,1),'n') then do
  54.     'Message "Invalid cell reference"'
  55.     'DEFPUBSCREEN("Workbench")'
  56.     exit
  57. end
  58. /* Suppress Screen Redraw to Speed Things Up */
  59. 'Refresh 0'
  60.  
  61. /* Open a small output window on tcalc screen*/
  62. fo=0
  63. CR='0a'x
  64. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  65. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  66.      call writeln(6Info, DisplayMsg)
  67.     fo=1
  68. end
  69. else do
  70.     'Message "TCALC Screen not available for Progress messages"'
  71. end
  72. CALL DELAY(150)
  73.  
  74. /* Get cell references for top cell in each column */
  75. 'SelectCell' start_cell
  76. do col=start_col to end_col
  77.     'GetCursorPos'
  78.     top_cell.col=result
  79.     'Column 1'
  80. end
  81.  
  82. /* Get labels for later use on output */
  83. 'SelectCell' start_cell
  84. 'GetValue'
  85. testlabel=result
  86. testlabel=strip(testlabel)
  87. if datatype(testlabel,'n')=1 then do
  88.     labelflag=0
  89.     do x=1 to NCols
  90.     title.x="Column "||x
  91.     end
  92. end
  93. else do
  94.     labelflag=1
  95.     NRows=NRows-1
  96.     do x=1 to NCols
  97.     'GetValue'
  98.     str=result
  99.     title.x=translate(strip(str),"_"," ")
  100.     'Column 1'
  101.     end
  102. end
  103. if fo then call writech(6Info,"Progress...10 ")
  104.  
  105. /* Get data from cell range */
  106. col=start_col
  107. lav=0
  108. tot=0
  109. count.=0
  110. total.=0
  111. do x=1 to NCols
  112.     'SelectCell' top_cell.col
  113.     if labelflag=1 then 'CursorDown 1'
  114.     do y=1 to NRows
  115.         'GetValue'
  116.         valtest=result
  117.         valtest=strip(valtest)
  118.         if datatype(valtest,'n')=1 then do
  119.             'GetValue'
  120.             val=result
  121.             data.x.y=val
  122.             tot=tot+val
  123.             total.x=tot
  124.             count.x=1+count.x
  125.         end
  126.         'CursorDown 1'
  127.     end
  128.     col=col+1
  129.     tot=0
  130.     lav=0
  131.     val=0
  132. end
  133. if fo then call writech(6Info,"20 ")
  134.  
  135. /* Calculate Means */
  136. mean.=0
  137. do x=1 to NCols
  138.     mean.x=total.x/count.x
  139. end
  140. if fo then call writech(6Info,"30 ")
  141.  
  142. /* Calculate Sum of Squared Values */
  143. sumsq.=0
  144. do x =1 to NCols
  145.     do y=1 to count.x
  146.         sumsq.x=(sumsq.x)+(data.x.y)**2
  147.     end
  148. end
  149. if fo then call writech(6Info,"40 ")
  150.  
  151. /* Calculate Unbiased Estimate of Variance */
  152. Var=(((sumsq.1)-((total.1)**2/count.1))+((sumsq.2)-((total.2)**2/count.2)))/((count.1+count.2)-2)
  153. StDev=sqrt(Var)
  154. StErr=sqrt((VAR/count.1)+(VAR/count.2))
  155.  
  156. /* Calculate t ratio */
  157. t=((mean.1)-(mean.2))/sqrt((Var/count.1)+(Var/count.2))
  158. if fo then call writech(6Info,"50 ")
  159.  
  160. /* Calculate degrees of freedom */
  161. df=(count.1)+(count.2)-2
  162.  
  163. /* Calculate T distribution function */
  164.  
  165. P=PROBT(t,df)
  166. if t>=0 then P=1-P
  167. PT=P*2
  168. if fo then call writech(6Info,"60 ")
  169. /* Calculate T Critical */
  170.  
  171. TCRIT1=TCRITICAL(.95,df)
  172. if fo then call writech(6Info,"70 ")
  173. TCRIT2=TCRITICAL(.99,df)
  174. if fo then call writech(6Info,"80 ")
  175. TCRIT3=TCRITICAL(.975,df)
  176. if fo then call writech(6Info,"90 ")
  177. TCRIT4=TCRITICAL(.995,df)
  178. if fo then do
  179.     call writeln(6Info,"100")
  180.     call writeln(6Info,"Writing output to window...")
  181. end
  182. /* Output */
  183. 'SelectCell' out_cell
  184. 'ColumnWidth 25' 
  185. 'Put "T-Test: Two Means for Independent Samples"'
  186. 'Column 1'
  187. 'CursorDown 1'
  188. do x=1 to NCols
  189.     'GetCursorPos'
  190.     first_cell.x=result
  191.     title=""""||title.x||""""
  192.     'Alignment 2'
  193.     'Put' title
  194.     'Column 1'
  195. end
  196. 'SelectCell' out_cell
  197. 'CursorDown 2'
  198. 'Put' "Count:"
  199. 'CursorDown 1'
  200. 'Put' "Sum:"
  201. 'CursorDown 1'
  202. 'Put' "Mean:"
  203. 'CursorDown 2'
  204. 'GetCursorPos'
  205. restart_cell=result
  206. 'Put "Pooled Variance:"'
  207. 'CursorDown 1'
  208. 'Put "Std. Deviation:"'
  209. 'cursorDown 1'
  210. 'Put "St. Error:"'
  211. 'CursorDown 1'
  212. 'Put "t:"'
  213. 'CursorDown 1'
  214. 'Put "d.f.:"'
  215. 'CursorDown 1'
  216. 'Put "P(T<=t) one-tail:"'
  217. 'CursorDown 1'
  218. 'Put "T-Critical (95%):"'
  219. 'CursorDown 1'
  220. 'Put "T-Critical (99%):"'
  221. 'CursorDown 1'
  222. 'Put "P(T<=t) two-tail:"'
  223. 'CursorDown 1'
  224. 'Put "T-Critical (95%):"'
  225. 'CursorDown 1'
  226. 'Put "T-Critical (99%):"'
  227. do x=1 to NCols
  228.     'SelectCell' first_cell.x
  229.     'CursorDown 1'
  230.     'Put' count.x
  231.     'CursorDown 1'
  232.     'Put' total.x
  233.     'CursorDown 1'
  234.     'Put' format(mean.x,,4)
  235.     end
  236. 'SelectCell' restart_cell
  237. 'Column 1'
  238. 'Put' format(Var,,4)
  239. 'CursorDown 1'
  240. 'Put' format(StDev,,4)
  241. 'CursorDown 1'
  242. 'Put' format(StErr,,4)
  243. 'CursorDown 1'
  244. 'Put' format(t,,4)
  245. 'CursorDown 1'
  246. 'Put' df
  247. 'CursorDown 1'
  248. 'Put' format(P,,6)
  249. 'CursorDown 1'
  250. 'Put' format(TCRIT1,,4)
  251. 'CursorDown 1'
  252. 'Put' format(TCRIT2,,4)
  253. 'CursorDown 1'
  254. 'Put' format(PT,,6)
  255. 'CursorDown 1'
  256. 'Put' format(TCRIT3,,4)
  257. 'CursorDown 1'
  258. 'Put' format(TCRIT4,,4)
  259. 'Refresh 1'
  260. 'Refresh 2'
  261. /*'Message' "Finished"*/
  262. /*indicate the main script is finished*/
  263. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  264. result=writeln(6Info, DisplayMsg)
  265. if result~=0 then do
  266.     /*Wait 3 seconds*/
  267.     CALL DELAY(150)
  268.     /* close window*/
  269.     result=close(6Info)
  270. end
  271. 'DEFPUBSCREEN("Workbench")'
  272. exit
  273.  
  274. /* Procedures */
  275.  
  276. cellrow: procedure
  277. do
  278.     parse arg cell
  279.     do charpos=2 to length(cell)
  280.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  281.     end
  282.     return 0
  283. end
  284. Return
  285.  
  286. cellcol: procedure
  287. do
  288.     parse arg cell
  289.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  290.     cell=upper(cell)
  291.     len=length(cell)
  292.     val=0
  293. do charpos=1 to len
  294.     if datatype(substr(cell,charpos,1),n) then
  295.     do cell=reverse(substr(cell,1,charpos-1))
  296.     do x=1 to length(cell)
  297.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  298.     end
  299.     return val
  300.     end
  301.     end
  302.     return 0
  303. end
  304. Return
  305. /* It is important to put the exposed array at the end of the next line */
  306. Sort: procedure expose NCols NRows data.
  307. do x=1 to NCols
  308. L=(xtoy(2,int(log(NRows)/log(2))))-1
  309.     Do Until L<1
  310.     L=trunc(int(L/2))
  311.     Do J=1 to L
  312.             Do K=J+L To NRows By L
  313.             I=K
  314.             dumdat=data.x.I
  315.             Do while I>L
  316.                 y=I-L
  317.                 If data.x.y ~> dumdat then Leave
  318.                 data.x.I=data.x.y
  319.                 I=I-L
  320.             End
  321.             data.x.I=dumdat
  322.             End
  323.         End
  324.     End
  325. End
  326. Return
  327.  
  328. syntax:
  329.      if arg(1)='FAIL' then do
  330.         'Message "Library is unavailable."'
  331.         'DEFPUBSCREEN("Workbench")'
  332.         exit
  333.         end
  334.     'DEFPUBSCREEN("Workbench")'
  335.     exit
  336.  
  337. PROBT: PROCEDURE
  338.  
  339.     PARSE ARG TA,K1
  340.     A=.36338023
  341.     W=ATAN(TA/SQRT(K1))
  342.     S=SIN(W)
  343.     C=COS(W)
  344.     L=K1-2*INT(K1/2)
  345.     IF L=0 THEN DO
  346.         T1=S
  347.         IF K1=2 THEN DO
  348.         PO=.5*(1+T1)
  349.         RETURN PO
  350.         END
  351.         T2=S
  352.         J1=-1
  353.         J2=0
  354.         K2=(K1-2)/2
  355.     END
  356.     ELSE DO
  357.         T1=W
  358.         IF K1=1 THEN DO
  359.             T1=T1*(1-A*L)
  360.             PO=.5*(1+T1)
  361.             RETURN PO
  362.         END
  363.         T2=S*C
  364.         T1=T1+T2
  365.         IF K1=3 THEN DO
  366.             T1=T1*(1-A*L)
  367.             PO=.5*(1+T1)
  368.             RETURN PO
  369.         END
  370.         J1=0
  371.         J2=1
  372.         K2=(K1-3)/2
  373.     END
  374.     DO I=1 TO K2
  375.         J1=J1+2
  376.         J2=J2+2
  377.         T2=T2*C*C*J1/J2
  378.         T1=T1+T2
  379.     END
  380.     T1=T1*(1-A*L)
  381.     PO=.5*(1+T1)
  382. RETURN PO
  383.  
  384. TCRITICAL: PROCEDURE
  385.  
  386.     PARSE ARG PO,K1
  387.     AO=.0001
  388.     E=.005
  389.     E2=E+E
  390.     A=2*PO-1
  391.     IF K1=1 THEN DO
  392.         TO=TAN(1.57079633*A)
  393.         RETURN TO
  394.     END
  395.     IF K1=2 THEN DO
  396.         SN=SIGN(2/(1-A*A))
  397.         IF SN=-1 THEN TO=-1*(A*SQRT(ABS(2/(1-A*A))))
  398.         ELSE TO=A*SQRT(2/(1-A*A))
  399.         RETURN TO
  400.     END
  401.     P1=PO
  402.     Z=NORMALPP(PO)
  403.     W=Z*(1+2/(1+8*K1))
  404.     T3=K1*(EXP(W*W/K1)-1)
  405.     SELECT
  406.         WHEN Z<0 THEN TT=-1
  407.         WHEN Z=0 THEN TT=0
  408.     OTHERWISE TT=1
  409.     END
  410.     SN=SIGN(T3)
  411.     IF SN=-1 THEN T3=-1*(TT*SQRT(ABS(T3)))
  412.     ELSE T3=TT*SQRT(T3)
  413.     /*
  414.     LABELA:
  415.     TO=T3
  416.     PO=PROBT(TO,K1)
  417.     F1=PO-P1
  418.     TO=T3+E
  419.     PO=PROBT(TO,K1)
  420.     F2=PO
  421.     TO=T3-E
  422.     PO=PROBT(TO,K1)
  423.     F2=F2-PO
  424.     F2=F2/E2
  425.     T4=T3-F1/F2
  426.     IF ABS(T4-T3)>AO THEN DO
  427.         T3=T4
  428.         SIGNAL 'LABELA'
  429.     END
  430.     */
  431.     T4=0
  432.     DO FOREVER
  433.         TO=T3
  434.         PO=PROBT(TO,K1)
  435.         F1=PO-P1
  436.         TO=T3+E
  437.         PO=PROBT(TO,K1)
  438.         F2=PO
  439.         TO=T3-E
  440.         PO=PROBT(TO,K1)
  441.         F2=F2-PO
  442.         F2=F2/E2
  443.         T4=T3-F1/F2
  444.         IF ABS(T4-T3)>AO THEN T3=T4
  445.         ELSE LEAVE
  446.     END
  447.     TO=T4
  448.     
  449. RETURN TO
  450.  
  451. LOGGAMMA: PROCEDURE
  452.  
  453.     ARG XA
  454.     C1=76.18009173
  455.     C2=-86.50532033
  456.     C3=24.01409822
  457.     C4=-1.231739516
  458.     C5=.001208580
  459.     C6=-.000005364
  460.     C7=2.506628275
  461.     X1=XA-1
  462.     W=X1+5.5
  463.     W=(X1+.5)*LN(W)-W
  464.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  465.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  466.     L=W+LN(C7*S)
  467. RETURN L
  468.  
  469. NORMALPP: PROCEDURE
  470.  
  471.     ARG P0
  472.     A1=2.515517
  473.     A2=.802853
  474.     A3=.010328
  475.     A4=1.432788
  476.     A5=.189269
  477.     A6=.001308
  478.     Q=.5-ABS(P0-.5)
  479.     SN=SIGN(-2*LN(Q))
  480.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  481.     ELSE W=SQRT(-2*LN(Q))
  482.     W1=A1+W*(A2+A3*W)
  483.     W2=1+W*(A4+W*(A5+A6*W))
  484.     ZZ=W-W1/W2
  485.     SELECT
  486.         WHEN (P0-.5)<0 THEN TT=-1
  487.         WHEN (P0-.5)=0 THEN TT=0
  488.         otherwise TT=1
  489.     END
  490.     ZZ=ZZ*TT
  491. RETURN ZZ
  492.  
  493. Format:  procedure
  494.  
  495.      arg number, before, after
  496.      CallLine = SIGL
  497.      if ~datatype(CallLine, 'N') then CallLine = '??'
  498.  
  499.      /* Make sure we have a number as first (required) argument    */
  500.      if ~datatype(number, 'N') then do
  501.         if number = '' then
  502.            fc = 17     /* Wrong number of arguments           */
  503.         else
  504.            fc = 47     /* Arithmetic conversion error             */
  505.         signal FormatSyntaxError
  506.      end
  507.      num = number + 0
  508.      if before = '' & after = '' then
  509.         return num
  510.      else do
  511.         parse var num integer '.' fraction
  512.         if before = '' then before = length(integer)
  513.         if after = '' then after = length(fraction)
  514.         if ~datatype(before, N) | ~datatype(after, N) then
  515.            do fc = 18
  516.            signal FormatSyntaxError
  517.        end
  518.         if before < length(integer) then do
  519.            fc = 18
  520.            signal FormatSyntaxError
  521.         end
  522.         if after ~= length(fraction) then do
  523.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  524.            integer = integer + (fraction % 1)
  525.            fraction = substr(fraction, 3)
  526.         end
  527.         if fraction >= 0 then
  528.            return right(integer, before)'.'fraction
  529.         else
  530.            return right(integer, before)
  531.      end
  532.  
  533.  FormatSyntaxError:
  534.         if show('F', STDERR) then
  535.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  536.         else
  537.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  538.         'Message' mess
  539.         parse source Func .
  540.         if Func = 'FUNCTION' then do
  541.         'DEFPUBSCREEN("Workbench")'
  542.            exit "Err"
  543.         end
  544.         else do
  545.         'DEFPUBSCREEN("Workbench")'
  546.            exit 10
  547.         end